module radar_netcdf_mod

  use container
  use flogger
  use string
  use data_translators_utils_mod
  use radar_mod
  use netcdf

  implicit none

  private

  public radar_netcdf_write

  interface
    subroutine add_string_vars(file_path, time_reference, sweep_mode, follow_mode, prt_mode) bind(c, name='add_string_vars')
      use iso_c_binding
      character(c_char), intent(in) :: file_path
      character(c_char), intent(in) :: time_reference
      character(c_char), intent(in) :: sweep_mode
      character(c_char), intent(in) :: follow_mode
      character(c_char), intent(in) :: prt_mode
    end subroutine add_string_vars
  end interface

contains

  subroutine radar_netcdf_write(file_path, stations, records)

    character(*), intent(inout) :: file_path
    type(hash_table_type), intent(in) :: stations
    type(linked_list_type), intent(in) :: records

    integer, parameter :: max_record_per_station = 10
    integer elevation_size, azimuth_size, bin_size

    type(linked_list_iterator_type) record_it
    character(30) units
    character(100) standard_name, long_name
    integer ncid, ierr, iscan, iradial, imoment, ibin
    integer max_range_size
    integer sweep_dimid
    integer time_coverage_start_varid
    integer time_coverage_end_varid
    integer lat_varid
    integer lon_varid
    integer alt_varid
    integer altitude_agl_varid
    integer sweep_grpid
    integer time_dimid
    integer range_dimid
    integer freq_dimid
    integer sweep_number_varid
    integer time_varid
    integer range_varid
    integer freq_varid
    integer ray_index_varid
    integer sweep_fixed_angle_varid
    integer azimuth_varid
    integer elevation_varid
    integer nyq_speed_varid
    integer field_varid

    if (file_path == '') file_path ='radar.nc'

    record_it = linked_list_iterator(records)
    do while (.not. record_it%ended())
      select type (record => record_it%value)
      type is (radar_record_type)
        ierr = nf90_create(file_path, nf90_netcdf4, ncid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'Conventions', 'Cf/Radial')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'version', '2.0')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'title', 'Radar base data')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'institution', 'CMA MOC')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'references', 'N/A')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'source', 'N/A')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'history', 'N/A')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'comment', 'Converted by data-translators')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'instrument_name', 'X-band radar')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'site_name', record%station%id)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, nf90_global, 'scan_name', record%scan_type)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_def_dim(ncid, 'sweep', size(record%scans), sweep_dimid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_def_var(ncid, 'time_coverage_start', nf90_double, varid=time_coverage_start_varid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, time_coverage_start_varid, 'units', 'seconds since 1970-01-01T00:00:00')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_def_var(ncid, 'time_coverage_end', nf90_double, varid=time_coverage_end_varid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, time_coverage_end_varid, 'units', 'seconds since 1970-01-01T00:00:00')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_def_var(ncid, 'latitude', nf90_double, varid=lat_varid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, lat_varid, 'units', 'degrees_north')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_def_var(ncid, 'longitude', nf90_double, varid=lon_varid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, lon_varid, 'units', 'degrees_east')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_def_var(ncid, 'altitude', nf90_double, varid=alt_varid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, alt_varid, 'units', 'meters')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_def_var(ncid, 'altitude_agl', nf90_double, varid=altitude_agl_varid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_att(ncid, altitude_agl_varid, 'units', 'meters')
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_def_var(ncid, 'sweep_fixed_angle', nf90_float, sweep_dimid, sweep_fixed_angle_varid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ! Find out maximum range size.
        max_range_size = 0
        do iscan = 1, size(record%scans)
          max_range_size = max(max_range_size, size(record%scans(iscan)%radials(1)%moments(1)%bins))
        end do

        ! Sweep groups
        do iscan = 1, size(record%scans)
          associate (scan => record%scans(iscan))
          ierr = nf90_def_grp(ncid, 'sweep_' // to_str(iscan, pad_zeros=2), sweep_grpid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_def_var(sweep_grpid, 'sweep_number', nf90_int, varid=sweep_number_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ! Define time dimension and variable.
          ierr = nf90_def_dim(sweep_grpid, 'time', scan%nradial, time_dimid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_def_var(sweep_grpid, 'time', nf90_double, time_dimid, time_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, time_varid, 'standard_name', 'time')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, time_varid, 'units', 'seconds since 1970-01-01T00:00:00')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ! Define range dimension and variable.
          ierr = nf90_def_dim(sweep_grpid, 'range', max_range_size, range_dimid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_def_var(sweep_grpid, 'range', nf90_double, range_dimid, range_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, range_varid, 'standard_name', 'projection_range_coordinate')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, range_varid, 'long_name', 'range to measurement volume')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, range_varid, 'units', 'meters')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, range_varid, 'spacing_is_constant', 'true')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, range_varid, 'meters_to_center_of_first_gate', scan%log_resolution*0.5)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, range_varid, 'meters_between_gates', scan%log_resolution)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, range_varid, 'axis', 'radial_range_coordinate')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ! Define frequency variable.
          ierr = nf90_def_dim(sweep_grpid, 'frequency', 1, freq_dimid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_def_var(sweep_grpid, 'frequency', nf90_real, freq_dimid, freq_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, freq_varid, 'units', 's-1')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ! Define ray_index variable.
          ierr = nf90_def_var(sweep_grpid, 'ray_index', nf90_int, time_dimid, ray_index_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ! Define sweep_fixed_angle variable.
          ierr = nf90_def_var(sweep_grpid, 'sweep_fixed_angle', nf90_real, varid=sweep_fixed_angle_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, sweep_fixed_angle_varid, 'units', 'degrees')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ! Define azimuth variable.
          ierr = nf90_def_var(sweep_grpid, 'azimuth', nf90_real, time_dimid, azimuth_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, azimuth_varid, 'standard_name', 'ray_azimuth_angle')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, azimuth_varid, 'long_name', 'azimuth angle from true north')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, azimuth_varid, 'units', 'degrees')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, azimuth_varid, 'axis', 'radial_azimuth_coordinate')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ! Define elevation variable.
          ierr = nf90_def_var(sweep_grpid, 'elevation', nf90_real, time_dimid, elevation_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, elevation_varid, 'standard_name', 'ray_elevation_angle')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, elevation_varid, 'long_name', 'elevation angle from horizontal plane')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, elevation_varid, 'units', 'degrees')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, elevation_varid, 'axis', 'radial_elevation_coordinate')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ! Define nyquist_velocity variable.
          ierr = nf90_def_var(sweep_grpid, 'nyquist_velocity', nf90_real, time_dimid, nyq_speed_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_att(sweep_grpid, nyq_speed_varid, 'units', 'meters per second')
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          do imoment = 1, size(scan%radials(1)%moments)
            associate (moment => scan%radials(1)%moments(imoment))
            select case (moment%name)
            case ('dBZ')
              units = 'dBZ'
              standard_name = 'equivalent_reflectivity_factor'
              long_name = 'reflectivity'
            case ('V')
              units = 'meters per second'
              standard_name = 'radial_velocity_of_scatters_away_from_instrument'
              long_name = 'mean doppler velocity'
            case default
              cycle
            end select
            ierr = nf90_def_var(sweep_grpid, moment%name, nf90_short, [range_dimid,time_dimid], field_varid)
            call handle_netcdf_error(ierr, __FILE__, __LINE__)

            ierr = nf90_put_att(sweep_grpid, field_varid, 'standard_name', standard_name)
            call handle_netcdf_error(ierr, __FILE__, __LINE__)

            ierr = nf90_put_att(sweep_grpid, field_varid, 'long_name', long_name)
            call handle_netcdf_error(ierr, __FILE__, __LINE__)

            ierr = nf90_put_att(sweep_grpid, field_varid, 'units', units)
            call handle_netcdf_error(ierr, __FILE__, __LINE__)

            ierr = nf90_put_att(sweep_grpid, field_varid, 'scale_factor', moment%scale_factor)
            call handle_netcdf_error(ierr, __FILE__, __LINE__)

            ierr = nf90_put_att(sweep_grpid, field_varid, 'add_offset', moment%add_offset)
            call handle_netcdf_error(ierr, __FILE__, __LINE__)
            end associate
          end do
          end associate
        end do

        ierr = nf90_enddef(ncid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ! Write root variables.
        ierr = nf90_put_var(ncid, time_coverage_start_varid, record%scans(1)%radials(1)%time%timestamp())
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_var(ncid, time_coverage_end_varid, record%scans(size(record%scans))%radials(size(record%scans(size(record%scans))%radials))%time%timestamp())
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_var(ncid, lat_varid, record%station%lat)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_var(ncid, lon_varid, record%station%lon)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_var(ncid, alt_varid, record%station%z + record%station%antenna_hgt)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_var(ncid, altitude_agl_varid, record%station%antenna_hgt)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_inq_varid(ncid, 'sweep_fixed_angle', sweep_fixed_angle_varid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ierr = nf90_put_var(ncid, sweep_fixed_angle_varid, record%scans(:)%elevation)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        ! Write sweep variables in each group.
        do iscan = 1, size(record%scans)
          associate (scan => record%scans(iscan))
          ierr = nf90_inq_ncid(ncid, 'sweep_' // to_str(iscan, pad_zeros=2), sweep_grpid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_inq_varid(sweep_grpid, 'sweep_number', sweep_number_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_var(sweep_grpid, sweep_number_varid, iscan-1)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_inq_varid(sweep_grpid, 'frequency', freq_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_var(sweep_grpid, freq_varid, record%station%freq)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_inq_varid(sweep_grpid, 'sweep_fixed_angle', sweep_fixed_angle_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_put_var(sweep_grpid, sweep_fixed_angle_varid, scan%elevation)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_inq_varid(sweep_grpid, 'time', time_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_inq_varid(sweep_grpid, 'range', range_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_inq_varid(sweep_grpid, 'azimuth', azimuth_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_inq_varid(sweep_grpid, 'elevation', elevation_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_inq_varid(sweep_grpid, 'nyquist_velocity', nyq_speed_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          ierr = nf90_inq_varid(sweep_grpid, 'ray_index', ray_index_varid)
          call handle_netcdf_error(ierr, __FILE__, __LINE__)

          do iradial = 1, scan%nradial
            associate (radial => scan%radials(iradial))
            ierr = nf90_put_var(sweep_grpid, ray_index_varid, [iradial], start=[iradial], count=[1])
            call handle_netcdf_error(ierr, __FILE__, __LINE__)

            ierr = nf90_put_var(sweep_grpid, time_varid, [radial%time%timestamp()], start=[iradial], count=[1])
            call handle_netcdf_error(ierr, __FILE__, __LINE__)

            ierr = nf90_put_var(sweep_grpid, azimuth_varid, [radial%azimuth], start=[iradial], count=[1])
            call handle_netcdf_error(ierr, __FILE__, __LINE__)

            ierr = nf90_put_var(sweep_grpid, elevation_varid, [radial%elevation], start=[iradial], count=[1])
            call handle_netcdf_error(ierr, __FILE__, __LINE__)

            ierr = nf90_put_var(sweep_grpid, nyq_speed_varid, [scan%nyq_speed], start=[iradial], count=[1])
            call handle_netcdf_error(ierr, __FILE__, __LINE__)
            end associate
          end do
          ! Write range variable.
          do ibin = 1, size(scan%radials(1)%moments(1)%bins)
            ierr = nf90_put_var(sweep_grpid, range_varid, [scan%log_resolution*(ibin-0.5)], start=[ibin], count=[1])
            call handle_netcdf_error(ierr, __FILE__, __LINE__)
          end do
          ! Write field variables.
          do iradial = 1, scan%nradial
            associate (radial => scan%radials(iradial))
            do imoment = 1, size(radial%moments)
              associate (moment => radial%moments(imoment))
              if (moment%name /= 'dBZ' .and. moment%name /= 'V') cycle
              ierr = nf90_inq_varid(sweep_grpid, moment%name, field_varid)
              call handle_netcdf_error(ierr, __FILE__, __LINE__)

              ierr = nf90_put_var(sweep_grpid, field_varid, (moment%bins - moment%add_offset) / moment%scale_factor, start=[1,iradial], count=[size(moment%bins),1])
              call handle_netcdf_error(ierr, __FILE__, __LINE__)
              end associate
            end do
            end associate
          end do
          end associate
        end do

        ierr = nf90_close(ncid)
        call handle_netcdf_error(ierr, __FILE__, __LINE__)

        call add_string_vars(trim(file_path) // char(0), &
          record%time%isoformat() // char(0), &
          'elevation_surveillance' // char(0), &
          'none' // char(0), &
          'fixed' // char(0))
      end select
      call record_it%next()
    end do

  end subroutine radar_netcdf_write

end module radar_netcdf_mod
